class verletBrownian
{
private:
	caja *box;
	///Paso de tiempo Dt
	double Dt;
	///Dt2
	double Dt2;
	///Dt/2
	double Dth;
	///Dt2/0.5
	double Dt2h;
	///Temperatura
	double temp;
	///Coeficiente de difusion (gamma=kBT/m*Diff)
	double diff;
	#ifndef POLYDISPERSE
    ///Coeficiente de friccion
	double gamma;
	///xsd = gamma*dt
	double xsd;
    ///c0 = exp[-gamma*dt]
    double c0;
    ///c1=(1-c0)/(gamma*dt)
    double c1;
    ///c2=(1-c1)/(gamma*dt)
    double c2;
    ///c1dt = c1*dt, c2dt=c2*dt
    double c1dt,c2dt;
    double c2dt2;
    ///c1mc2dt=(c1-c2)*dt
    double c1mc2dt;
    ///Para generar los desplazamientos y velocidades aleatorias
    double sigmaR, sigmaV;
    double corrrV,umcorrV2SQ;
    #endif

    #ifdef POLYDISPERSE
	double *gamma;
	///xsd = gamma*dt
	double *xsd;
    ///c0 = exp[-gamma*dt]
    double *c0;
    ///c1=(1-c0)/(gamma*dt)
    double *c1;
    ///c2=(1-c1)/(gamma*dt)
    double *c2;
    ///c1dt = c1*dt, c2dt=c2*dt
    double *c1dt,*c2dt;
    double *c2dt2;
    ///c1mc2dt=(c1-c2)*dt
    double *c1mc2dt;
    ///Para generar los desplazamientos y velocidades aleatorias
    double *sigmaR, *sigmaV;
    double *corrrV, *umcorrV2SQ;
    #endif

    double **drg, **dvg;

public:

    verletBrownian(caja *bx, double xDt, double xtemp, double xdiff)
    {
        box=bx;
        temp=xtemp;
        diff=xdiff;
        Dt=xDt;

        #ifndef POLYDISPERSE
        //gamma = kB*temp/diff;
        gamma = temp/diff;
        xsd=gamma*Dt;  //(kBT·t/(m·Dif))
        c0=exp(-xsd);
        c1=(1.0-c0)/xsd;
        c2=(1.0-c1)/xsd;
        c1dt=c1*Dt;
        c2dt=c2*Dt;
        c2dt2=c2*Dt*Dt;
        c1mc2dt=(c1-c2)*Dt;
        sigmaR= sqrt(Dt*Dt*kB*temp/xsd*(2.0-(3.0-4.0*exp(-xsd)+exp(-2.0*xsd))/xsd));
        sigmaV= sqrt(kB*temp*(1.0-exp(-2.0*xsd)));
        corrrV=Dt*kB*temp/xsd*(1.0-exp(-xsd))*(1.0-exp(-xsd))/(sigmaR*sigmaV);
        umcorrV2SQ=sqrt(1.0-corrrV*corrrV);
        #endif

        #ifdef POLYDISPERSE
        gamma = new double[bx->npart];
        xsd = new double[bx->npart];
        c0=new double[bx->npart];
        c1=new double[bx->npart];
        c2=new double[bx->npart];
        c1dt=new double[bx->npart];
        c2dt=new double[bx->npart];
        c2dt2=new double[bx->npart];
        c1mc2dt=new double[bx->npart];
        sigmaR=new double[bx->npart];
        sigmaV=new double[bx->npart];
        corrrV=new double[bx->npart];
        umcorrV2SQ=new double[bx->npart];

        for (int i=0;i<bx->npart;i++)
        {
            gamma[i]=kB*temp/diff/bx->part[i].mass;
            xsd[i]=gamma[i]*Dt;
            c0[i]=exp(-xsd[i]);
            c1[i]=(1.0-c0[i])/xsd[i];
            c2[i]=(1.0-c1[i])/xsd[i];
            c1dt[i]=c1[i]*Dt;
            c2dt[i]=c2[i]*Dt;
            c2dt2[i]=c2[i]*Dt*Dt;
            c1mc2dt[i]=(c1[i]-c2[i])*Dt;
            sigmaR[i]= sqrt(Dt*Dt*kB*temp/bx->part[i].mass/xsd[i]*(2.0-(3.0-4.0*exp(-xsd[i])+exp(-2.0*xsd[i]))/xsd[i]));
            sigmaV[i]= sqrt(kB*temp/bx->part[i].mass*(1.0-exp(-2.0*xsd[i])));
            corrrV[i]=Dt*kB*temp/bx->part[i].mass/xsd[i]*(1.0-exp(-xsd[i]))*(1.0-exp(-xsd[i]))/(sigmaR[i]*sigmaV[i]);
            umcorrV2SQ[i]=sqrt(1.0-corrrV[i]*corrrV[i]);
        }

        #endif
        drg = new double*[box->npart];
        dvg= new double *[box->npart];
        for (int i=0;i<box->npart;i++)
        {
            drg[i]=new double[3];
            dvg[i]=new double[3];
        }
    }

    ~verletBrownian()
    {
        for (int i=0;i<box->npart;i++)
        {
            delete drg[i];
            delete dvg[i];
        }
        delete []drg;
        delete []dvg;

        #ifdef POLYDISPERSE
        delete []gamma;
        delete []xsd;
        delete []c0;
        delete []c1;
        delete []c2;
        delete []c1dt;
        delete []c2dt;
        delete []c2dt2;
        delete []c1mc2dt;
        delete []sigmaR;
        delete []sigmaV;
        delete []corrrV;
        delete []umcorrV2SQ;
        #endif
    }
	///Motor Verlet-Langevin
	void  verletLangevin(void);
	///Generador variables aleatorias gaussianas
	double gaussianRand(void);
	///Componentes aleatorias posicion-velocidad
	void gaussian(void);


};

void verletBrownian::gaussian(void)
{
    double xs1, xs2;
    for(int i=0;i<box->npart;i++)
    {
        for(int j=0;j<3;j++)
        {
            xs1=gaussianRand();
            xs2=gaussianRand();
            #ifndef POLYDISPERSE
            ///variable posicion aleatoria (ver Allen- Tildeslay)
            drg[i][j]=sigmaR*xs1;
            ///variable velocidad aleatoria (ver Allen- Tildeslay)
            dvg[i][j]=sigmaV*(corrrV*xs1+umcorrV2SQ*xs2);
            #endif
            #ifdef POLYDISPERSE
            ///variable posicion aleatoria (ver Allen- Tildeslay)
            drg[i][j]=sigmaR[i]*xs1;
            ///variable velocidad aleatoria (ver Allen- Tildeslay)
            dvg[i][j]=sigmaV[i]*(corrrV[i]*xs1+umcorrV2SQ[i]*xs2);
            #endif
        }
    }
}

double verletBrownian::gaussianRand(void)
{
     double xsi =0.0;
     for(int i=0;i<12;i++)
         xsi+=double(rand())/RAND_MAX;
     xsi-=6.0;
     return xsi;
}

void  verletBrownian::verletLangevin(void)
{
    double svxy, svxz, svyz, svzz;
    svxy=svxz=svyz=svzz=0.0;
    ///Creamos las componentes aleatorias
    gaussian();
    /// Actualizamos las posiciones de las particulas (ver Allen- Tildeslay)
    for(int i=0;i<box->npart;i++)
    {
        for (int j=0;j<3;j++)
        {
            #ifndef POLYDISPERSE
            box->part[i].pos[j]+=c1*box->part[i].vel[j]*Dt+
            c2*box->part[i].force[j]*Dt*Dt +drg[i][j];
            ///Primera actualizacion de las velocidades
            box->part[i].vel[j]=c0*box->part[i].vel[j]+(c1-c2)*box->part[i].force[j]*Dt+dvg[i][j];
            #endif

            #ifdef POLYDISPERSE
            box->part[i].pos[j]+=c1[i]*box->part[i].vel[j]*Dt+
                c2[i]*box->part[i].force[j]*Dt*Dt +drg[i][j];
            ///Primera actualizacion de las velocidades
            box->part[i].vel[j]=c0[i]*box->part[i].vel[j]+(c1[i]-c2[i])*box->part[i].force[j]*Dt+dvg[i][j];
            #endif
        }
        box->boundary(i);
    }

    #ifdef CONFINED
	interaccionConfined(box);
	#endif
	#ifndef CONFINED
	interaccionSimple(box);
	#endif
	double tauT=1.2;
    double alpha;
//    if(box->nacc==0)
//        alpha=1.0;
//    else
//        alpha=sqrt(1.0+Dt*(box->temp0/box->temp-1.0)/tauT) ;
//    alpha=1.0;
    //double alpha=sqrt(box->temp0/box->temp) ;
    //cout<<box->temp0<<"  "<<box->temp<<endl;
    ///Actualizamos las velocidades de las particulas (ver Allen- Tildeslay)
    double parke=0.0;
    for(int i=0;i<box->npart;i++)
    {
        for(int j=0;j<3;j++)
        {
            #ifndef POLYDISPERSE
            box->part[i].vel[j]+=c2*box->part[i].force[j]*Dt;
            parke+=box->part[i].vel[j]*box->part[i].vel[j];
            #endif
            #ifdef POLYDISPERSE
            box->part[i].vel[j]+=c2[i]*box->part[i].force[j]*Dt;
            parke+=box->part[i].mass*box->part[i].vel[j]*box->part[i].vel[j];
            #endif
            //box->part[i].vel[j]*=alpha;

        }

        #ifdef CALC_STRESS
			///Calculamos la componente de velocidad del tensor de esfuerzos
			#ifndef POLYDISPERSE
                svxy+=(box->part[i].vel[0])*(box->part[i].vel[1]);
                svxz+=(box->part[i].vel[0])*(box->part[i].vel[2]);
                svyz+=(box->part[i].vel[1])*(box->part[i].vel[2]);
                svzz+=(box->part[i].vel[2])*(box->part[i].vel[2]);
			#endif

			#ifdef POLYDISPERSE
                svxy+=(box->part[i].mass)*(box->part[i].vel[0])*(box->part[i].vel[1]);
                svxz+=(box->part[i].mass)*(box->part[i].vel[0])*(box->part[i].vel[2]);
                svyz+=(box->part[i].mass)*(box->part[i].vel[1])*(box->part[i].vel[2]);
                svzz+=(box->part[i].mass)*(box->part[i].vel[2])*(box->part[i].vel[2]);
			#endif
		#endif
    }

    #ifdef CALC_STRESS
		///Actualizamos el valor de la componente de velocidad del tensor de esfuerzos
		box->svxy=svxy/box->vol;
		box->svxz=svxz/box->vol;
		box->svyz=svyz/box->vol;
		box->svzz=svzz/box->vol;
	#endif



	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/double(box->npart);
	///Calculamos la temperatura de la caja
	box->temp=parke/3.0/double(box->npart);
	///Calculamos la presion de la caja
	box->presion=box->rho*box->temp+box->virial/box->vol;
}
